Ladislas Nalborczyk & Thierry Phénix
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 & 4: Modèle de régression linéaire
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Modèles multi-niveaux
Cours 9: Introduction à brms
Cours 10: Data Hackaton
Un ensemble de données / mesures observées à analyser :
\( d_1, d_2, \dots, d_N \)
Un modèle génératif défini par un ensemble de paramètres
\( \theta_1, \theta_2, \dots, \theta_K \)
Une connaissance a priori de ces paramètres
Nous avons les observations
The Bayesian approach places uncertainty not in the observations but rather in one’s lack of knowledge. For a Bayesian, the observed data are not uncertain—you observed what you observed.
Feinberg and Gonzalez (2012)
Un ensemble de données / mesures observées à analyser :
\( d_1, d_2, \dots, d_N \)
Un modèle génératif défini par un ensemble de paramètres
\( \theta_1, \theta_2, \dots, \theta_K \)
Une connaissance a priori de ces paramètres
Nous cherchons
Ces observations et ces paramètres sont regroupés dans un même modèle bayésien :
la distribution conjointe \( P(d_{1:N}~\theta_{1:K}) \).
\[ P(d_{1:N}~\theta_{1:K})~=~\color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})} \]
\[ P(d_{1:N}~\theta_{1:K})~=~\color{purple}{P(\theta_{1:K}~|~d_{1:N})}~\color{green}{P(d_{1:N})} \]
Ces deux décompositions sont égales !!!
On infére les nouvelles valeurs des paramètres à partir des données
\[ \color{purple}{P(\theta_{1:K}~|~d_{1:N})} = \frac{\color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})}}{\color{green}{P(d_{1:N})}} \]
Ou encore
\[ \color{purple}{P(\theta_{1:K}~|~d_{1:N})} \propto \color{orangered}{P(d_{1:N}~|~\theta_{1:K})}~\color{steelblue}{P(\theta_{1:K})} \]
Passer d'une connaissance a priori sur les paramètres \( \color{steelblue}{P(\theta_{1:K})} \)
À une connaissance à postériori \( \color{purple}{P(\theta_{1:K}~|~d_{1:N})} \)
1. Définir le modèle
2. Mettre à jour le modèle
3. Interpréter la distribution postérieure (répondre à la question posée)
4. Vérification prédictive du modèle
Illustrer les différents aspects de cette démarche
à l'aide d'un modèle simple
le modèle Beta-Binomial
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 & 4: Modèle de régression linéaire
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Modèles multi-niveaux
Cours 9: Introduction à brms
Cours 10: Data Hackaton
library(tidyverse)
set.seed(1234567)
sample(x = c(0, 1), size = 50, replace = TRUE, prob = c(0.4,0.6)) %>%
data.frame %>%
mutate(x = seq_along(.), y = cumsum(.) / seq_along(.) ) %>%
ggplot(aes(x = x, y = y), log = "y") + geom_line(lwd = 1, col = "steelBlue") +
geom_hline(yintercept = 0.5, lty = 3) +
xlab("nombre de lancers") + ylab("proportion de faces") +
ylim(0, 1) + theme_bw(base_size = 20)
Modèle d'un lancer unique
Modèle d'un lancer unique
Depuis Bernoulli, on sait calculer la probabilité du résultat d'un lancer, du moment que l'on connait le biais de la pièce
\[ P(Y=y~|~\theta)=\theta^y~(1-\theta)^{(1-y)} \]
avec \( \theta \) la probabilité de l'évènement “le résultat du lancer est \( \textit{Face} \)”
Modèle d'une suite de lancers (processus de Bernoulli)
Pour une série de lancers \( \{Y_i\} \) indépendants et identiquement distribués (i.i.d.)
\[ P(\{Y_i=y_i\}~|~\theta)=\theta^z~(1-\theta)^{(N-z)} \]
avec \( z=\sum_i y_i \) (nombre de succès) et N nombre d'observations
La fonction de likelihood s'écrit :
\[ L(\theta~|~z~N)=P(z~|~\theta~N)=\theta^z~(1-\theta)^{(N-z)} \]
c'est une fonction de la variable \( \theta \)
Exemple :
Représentation graphique de la fonction de vraissemblance de \( \theta \) pour 10 lancers et 7 Faces
NbrTrial = 100
N = 10
nbrSuccess = 7
theta = seq(from = 0, to = 1, by = 1/NbrTrial)
likelihood = (theta^(nbrSuccess)) * (1 - theta)^(N - nbrSuccess)
data.frame(theta, likelihood) %>%
ggplot(aes(x = theta, y = likelihood)) +
geom_area(color="orangered", fill = "orangered", alpha =.4) +
theme_bw(base_size = 20) + theme(legend.position="none",
text = element_text(size=30))
Aspect sémantique, rendre compte ou modéliser :
Aspect mathématique, pour une solution entièrement analytique :
Un candidat : La distribution beta
La distribution beta:
\[ \begin{align} p(\theta~|~a,b) &= beta(\theta~|~a,b) \\ &= \theta^{(a-1)}(1-\theta)^{(b-1)}/B(a,b) \\ &\propto \theta^{(a-1)}(1-\theta)^{(b-1)} \end{align} \]
où \( a \) et \( b \) sont deux paramètres tels que \( a\geq0~et~b\geq0 \) et \( B(a,b) \) est une constante de normalisation
Attention : Ne pas confondre dans R : la fonction de \( \theta \) et la constante (indépendante de \( \theta \))
a = 2
b = 3
theta = seq( from= 0 , to= 1 , length.out= 7 )
dbeta(theta,a,b)
[1] 0.0000000 1.3888889 1.7777778 1.5000000 0.8888889 0.2777778 0.0000000
beta(a,b)
[1] 0.08333333
Vérifions les contraintes mathématiques
Le résultat du produit d'une distribution beta et d'une fonction de vraissemblance de Bernoulli est proportionnel à une distribution beta.
On dit que la distribution beta est la distribution conjugée de la fonction de vraissemblance de Bernoulli.
Soit un prior définit par \( ~~~\color{steelblue}{p(\theta~|~a,b) \propto \theta^{(a-1)}(1-\theta)^{(b-1)}} \)
Soit une fonction de vraissemblance associée à z “Face” pour N lancers \( ~~\color{orangered}{P(\{y_i\}~|~\theta)=\theta^z~(1-\theta)^{(N-z)}} \)
\[ \begin{align} \color{purple}{p(\theta~|~\{y_i\})} &\propto \color{orangered}{P(\{y_i\}~|~\theta)}~\color{steelblue}{p(\theta~|~a,b)} &&\mbox{inférence Bayésienne} \\ &\propto \color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\theta^{(a-1)}(1-\theta)^{(b-1)}} &&\mbox{application des formules précédentes} \\ &\propto \color{purple}{\theta^{(z+a-1)}~(1-\theta)^{(N-z+b-1)}} &&\mbox{en regroupant les termes identiques} \\ &\propto \color{purple}{\theta^{(A-1)}(1-\theta)^{(B-1)}} &&\mbox{avec } A=z+a \mbox{ et } B=N-z+b \end{align} \]
Conclusions : la distribution a priori et la distribution a posteriori ont la même forme
Vérifions les aspects sémantiques
Vérifions les aspects sémantiques
Le niveau de certitude augmente avec la somme (\( \kappa = a+b \))
Le plus souvent, on dispose :
Le mode
Supposons que l'on ait une estimation de la valeur la plus probable \( \omega \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior
\[ \begin{align} a &= \omega(\kappa-2)+1 \\ b &= (1-\omega)(\kappa-2)+1 &&\mbox{pour } \kappa>2 \end{align} \]
Pour \( ~~\omega=0.65 \) et \( ~~\kappa=25 \) alors \( ~~~p(\theta)=beta(\theta~|~15.95, 9.05)~ \) (courbe orange)
Pour \( ~~\omega=0.65 \) et \( ~~\kappa=10 \) alors \( ~~~p(\theta)=beta(\theta~|~6.2, 3.8)~ \) (courbe bleue)
La moyenne
Supposons que l'on ait une estimation de la valeur moyenne \( \mu \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior
\[ \begin{align} a &= \mu \kappa \\ b &= (1-\mu) \kappa &&\mbox{ pour } \kappa>2 \end{align} \]
Pour \( \mu=0.8 \) et \( \kappa=7 \) alors \( p(\theta)=beta(\theta~|~5.6,1.4) \) (courbe orange)
Pour \( \mu=0.8 \) et \( \kappa=22 \) alors \( p(\theta)=beta(\theta~|~17.6,4.4) \) (courbe bleue)
La moyenne et la variance
Supposons que l'on ait une estimation a priori de la valeur moyenne \( \mu \) et de la variance \( \sigma \) du paramètre \( \theta \)
On peut définir les paramètres a et b de notre prior
\[ \begin{align} a &= \mu~(\mu~(1-\mu)~/~\sigma^2-1) \\ b &= (1-\mu)~\mu~(\mu~(1-\mu)~/~\sigma^2-1) ,~~~~~ \mbox{ pour } \sigma < 0.28867 \\ \end{align} \]
Pour \( ~~\mu=0.8 \) et \( ~~\sigma=0.15 \) alors \( ~~~p(\theta)=beta(\theta~|~5.92, 3.19) \) (orange)
Pour \( ~~\mu=0.8 \) et \( ~~\sigma=0.1 \) alors \( ~~~p(\theta)=beta(\theta~|~14.14, 7.61) \) (bleue)
Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :
\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]
Cas \( \mbox{ } N < a+b \)
Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :
\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]
Cas \( \mbox{ } N = a+b \)
Pour une distribution a priori \( \color{steelblue}{beta(\theta~|~a,b)} \), de moyenne \( \color{steelblue}{\frac{a}{a+b}} \), et un résultat de \( \color{orangered}{z} \) Face sur \( \color{orangered}{N} \) lancers, la moyenne de la distribution postérieure est :
\[ \frac{z+a}{N+a+b} = \color{orangered}{\frac{z}{N}}~\frac{N}{N+a+b}+\color{steelblue}{\frac{a}{a+b}}~\frac{a+b}{N+a+b} \]
Cas \( \mbox{ } N > a+b \)
There are two main desiderata for a mathematical description of data.
- First, the mathematical form should be comprehensible with meaningful parameters.
- The second desideratum for a mathematical description is that it should be descriptively adequate, which means, loosely, that the mathematical form should “look like” the data.
Kruschke (2015)
Cours 1 : Introduction à l'inférence bayésienne
Cours 2 : Modèle beta-binomial
Cours 3 & 4 : Modèle de régression linéaire
Cours 5 : Markov Chain Monte Carlo
Cours 6 : Modèle linéaire généralisé
Cours 7 : Comparaison de modèles
Cours 8 : Modèles multi-niveaux
Cours 9 : Introduction à brms
Cours 10 : Data Hackaton
The posterior distribution is always a compromise between the prior distribution and the likelihood function.
Kruschke (2015)
\[ Posterior = \frac{Likelihood~\times~Prior}{\color{orangered}{Marginal~Likelhood}} \]
le calcul de la “marginal likelihood”, \( \color{green}{~p(D)} \)
\[ \begin{align} \color{green}{p(D)} &= \int P(D,\theta) \, \mathrm d\theta &&\mbox{marginalisation sur le paramètre} \theta \\ \color{green}{p(D)} &= \color{green}{\int P(D~|~\theta)~P(\theta) \, \mathrm d\theta} && \mbox{Application de la règle du produit} \end{align} \]
Dans la littérature, on le trouve sous le nom de : “evidence”, “probabilité des données” ou “Marginal Likelihood”
\[ Posterior = \frac{Likelihood~\times~Prior}{\color{orangered}{Marginal~Likelhood}} \]
- Solution analytique \( ~~\longrightarrow~~ \) Utilisation d'une likelihood conjuguée du prior
- Solution discrètisée \( ~~\longrightarrow~~ \) Calcul de la solution sur un ensemble fini de points (Méthode Grid)
- Solution approchée \( ~~\longrightarrow~~ \) On échantillonne “intelligemment” sur l'espace conjoint des paramètres (Méthode MCMC)
Utilisation d'une likelihood conjuguée du prior : cas Beta-binomiale
\[
\begin{align}
\color{purple}{p(\theta~|~z,N)} &= \color{orangered}{P(z~|~\theta,N)}~\color{steelblue}{p(\theta)/\color{green}{P(z~|~N)}} &&\mbox{Équation de Bayes} \\
&= \color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\frac{\theta^{(a-1)}(1-\theta)^{(b-1)}}{B(a,b)}}/\color{green}{P(z~|~N)} &&\mbox{Application des formules précédetes} \\
&= \frac{\color{orangered}{\theta^z~(1-\theta)^{(N-z)}}~\color{steelblue}{\theta^{(a-1)}(1-\theta)^{(b-1)}}}{\color{steelblue}{B(a,b)}\color{green}{P(z~|~N)}} \\
&= \frac{\theta^{(\color{orangered}{z}+\color{steelblue}{a-1})}~(1-\theta)^{(\color{orangered}{N-z}+\color{steelblue}{b-1})}}{\color{steelblue}{B(a,b)}\color{green}{P(z~|~N)}} &&\mbox{En regroupant les termes identiques} \\
&= \frac{\theta^{(\color{orangered}{z}+\color{steelblue}{a-1})}~(1-\theta)^{(\color{orangered}{N-z}+\color{steelblue}{b-1})}}{B(\color{orangered}{z}+\color{steelblue}{a},\color{orangered}{N-z}+\color{steelblue}{b})} &&\mbox{Par identification à une fonction beta} \\
\end{align}
\]
L'avantage d'avoir un prior conjugué nous dispense du calcul du dénominateur
Distributions discrètes
Distributions continues
TROP CONTRAIGNANT !!!
Le modèle (likelihood + prior) doit être défini à partir de l'interprétation que l'on peut faire des paramètres et des distributions
et non pour faciliter les calculs
En affinant la grille on augmente le temps de calcul
Le “superordinateur” chinois Tianhe-2 réalise \( 33,8 x 10^{15} \) opérations par seconde
Si on considère qu'il réalise 1 opération par noeud de la grille, il lui faudrait \( 2.96 x 10^{13} \) pour parcourir la grille une fois !!!
Par comparaison, l'âge de l'univers est approximativement de \( (4,354 ± 0,012)×10^{17} \) secondes
trajLength = 250 # nbr of measure
theta = 1:7
ptheta = theta
trajectory = sample(theta, prob = ptheta, size = trajLength, replace = TRUE)
p_grid = seq( from=0 , to=1 , length.out=1000 )
a = b = 1
N = 9
z = 6
posterior = dbeta(p_grid, z+a-1, N-z+b-1)
p_grid = seq( from=0 , to=1 , length.out = 1000 )
prior = rep( 1 , 1000 )
likelihood = dbinom( 6 , size=8 , prob=p_grid )
posterior = likelihood * prior
posterior = posterior / sum(posterior)
samples = sample( p_grid , prob = posterior , size = 1e4 , replace = TRUE )
Méthode analytique
Méthode Grid
Méthodes par échantillonnage (MCMC)
La distribution postérieure est presque toujours représentée par un échantillon
Cours 1 : Introduction à l'inférence bayésienne
Cours 2 : Modèle beta-binomial
Cours 3 & 4 : Modèle de régression linéaire
Cours 5 : Markov Chain Monte Carlo
Cours 6 : Modèle linéaire généralisé
Cours 7 : Comparaison de modèles
Cours 8 : Modèles multi-niveaux
Cours 9 : Introduction à brms
Cours 10 : Data Hackaton
DESCRIPTION : A partir d'un échantillon d'une distribution postérieure, on peut calculer la moyenne, le mode et la médiane
library(rethinking)
chainmode( samples , adj=0.01 ) # pour le mode ou MAP
mean( samples ) # pour la moyenne
median( samples ) # pour la médiane
Exemple pour un prior uniforme, et 10 lancers successifs et 3 Faces
DESCRIPTION : A partir d'un échantillon d'une distribution postérieure, on peut calculer la moyenne, le mode et la médiane
library(rethinking)
chainmode( samples , adj=0.01 ) # pour le mode ou MAP
mean( samples ) # pour la moyenne
median( samples ) # pour la médiane
Exemple pour un prior uniforme, et 3 lancers successifs et 3 Faces
DÉCISION : Comment chosir la valeur qui résume le “mieux” la distribution postérieure
Utiliser des fonctions d'espérance de gain ou, ce qui revient au même, de perte attendue
# pour la moyenne
sapply( p_grid , function(d) sum( posterior*( d - p_grid )^2 ) )
# pour la médiane
sapply( p_grid , function(d) sum( posterior*abs( d - p_grid ) ) )
Définition: les valeurs du paramètre x contenues dans le HDI 89% sont telles que \( p(x)>W \) et
\[ \int_{x:p(x)>W} p(x) \, \mathrm dx=0.89 \]
En pratique - Package BEST
library(coda)
library(BEST)
set.seed(1234567)
iterMax = 10000
p_grid = seq( from=0 , to=1 , length.out=iterMax )
pTheta = dbeta( p_grid, 3, 10)
massVec = pTheta/sum(pTheta)
samples = sample( p_grid , size=1e4 , replace=TRUE , prob=pTheta )
#plotPost(samples, credMass = 0.85)
hdi( samples , credMass = 0.89 )
lower upper
0.05280528 0.39613961
attr(,"credMass")
[1] 0.89
En pratique - Package BEST
En pratique - Package rethinking
library(rethinking)
set.seed(1234567)
iterMax = 10000
p_grid = seq( from=0 , to=1 , length.out=iterMax )
pTheta = dbeta( p_grid, 3, 10)
massVec = pTheta / sum(pTheta)
samples = sample( p_grid , size = 1e4 , replace = TRUE , prob = pTheta )
# utilise le package CODA pour calculer HPDI
HPDI( samples, prob = 0.89)
|0.89 0.89|
0.05280528 0.39613961
On l'utilise pour “tester une hypothèse”
On considère une pièce. On la lance 200 fois et on obtient 115 Faces.
QUESTION : La pièce est-elle biaisée?
Nous construisons deux modèles et essayons de savoir lequel est le plus probable
\[ \begin{eqnarray*} \left\{ \begin{array}{ll} M_0~:~X \sim \mathcal{B}(N,0.5) & \quad \text{la pièce n'est pas biaisée}\\ M_1~:~X \sim \mathcal{B}(N,\theta) & \quad \text{la pièce est biaisée} \end{array}\right. \end{eqnarray*} \]
Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.
\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\color{orangered}{\frac{P(D~|~M_0)}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]
Dans cet exemple, le calcul de \( P(D~|~M_0) \) est facile
\( P(D~|~M_0) = C_{200}^{115}~(0.5)^{115}~(1-0.5)^{(200-115)} = 0.005955 \)
Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.
\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\frac{P(D~|~M_0)}{\color{orangered}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]
Pour calculer de \( P(D~|~M_1) \), c'est plus compliqué : \( \color{orangered}{~~~\textbf{On ne connait pas }\theta} \)
Le calcul de la vraissemblance dans ce cas se fait par marginalisation sur le paramètre \( \theta \) \[ \begin{align} P(D~|~M_1) &= \int_0^1 P(D~|~M_1, \theta) P(\theta) \, \mathrm d\theta \\ &= \int_0^1 C^z_N~\theta^z~(1-\theta)^{N-z}~\text{d}\theta \\ &= N+1 \end{align} \]
Le facteur de Bayes fait le rapport des vraissemblances des deux modèles.
\[ \frac{P(M_0~|~D)}{P(M_1~|~D)}=\color{orangered}{\frac{P(D~|~M_0)}{P(D~|~M_1)}}~\frac{P(M_0)}{P(M_1)} \]
Dans l'exemple,
\[ B_{01} = \frac{P(M_0~|~D)}{P(M_1~|~D)} = \frac{0.005955}{0.005} = 1.1971 \]
Le rapport de probabilités a augmenté de 20% après avoir pris connaissance des données.
La balance penche ici en faveur du modèle nul, \( B_{01} > 1 \)
BIC critère d’information bayésien
Schwarz (1977) propose une approximation du log de la vraissemblance:
La vraissemblance du modèle pas toujours calculable : \[ P(D~|~M)=\int_{\theta} P(D~|~M, \theta) \, \mathrm d\theta \]
Calcul d'une valeur approchée \[ ln(P(D~|~\theta)) \approx ln(L(\widehat{\theta})) - \frac{t}{2}~ln(N) \]
avec \( L(\widehat{\theta}) \) le maximum de vraissemblance
\( t \) le nombre de paramètres
\( N \) le nombre d'observations
Les deux rôles de la fonction de vraissemblance
On peut utiliser cette distribution de probabilité pour générer des données
Exemple
Générer 10000 valeurs à partir d'une loi binomiale basée sur 9 lancers et une probabilité de Face de 0,6 :
w <- rbinom(n = 10e4, size = 9, prob = 0.6)
Deux sources d'incertitude dans ces prédictions
Exemple
Générer 10000 valeurs à partir d'une loi binomiale basé sur 9 lancers et une probabilité de Face de 0,6
pour chaque valeur de p
w <- rbinom(n = 10e4, size = 9, prob = samples)